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Regions  Performance  of  Seismic-Acoustic  Sensor  Systems. 

The  author  thanks  Steven  Daly  and  Gary  DeCoff  of  CRREL  for  reviewing 
the  manuscript  of  this  report. 
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FORTRAN  SUBROUTINES  FOR  ZERO-PHASE  DIGITAL  FREQUENCY  FILTERS 

Donald  G.  Albert 


INTRODUCTION 

This  report  describes  FORTRAN  subroutines  that  can  be  used  to  design 
and  apply  zero-phase  recursive  frequency  filters  to  digital  data.  Four 
types  of  filters  are  included  and  are  described  bv  their  transfer  function 
characteristics  (Fig.  1):  low  pass  (LP),  high  pass  (HP),  bandpass  (BP)  and 
bandstop  (BS).  The  subroutines  allow  the  user  to  remove  noise  from  record¬ 
ed  data  or  to  enhance  a  selected  portion  of  the  signal.  The  performance  of 
the  filters  is  quite  good,  and  the  zero-nhase  property  of  the  subroutines 
that  applv  the  filters  ensures  that  no  unwanted  time  shifts  in  the  data 
will  be  produced. 

Digital  filtering  is  a  broad  field  with  a  vast  amount  of  literature, 
and  it  is  not  the  goal  of  this  report  to  describe  the  benefits  and  short¬ 
comings  of  various  filter  types  and  design  procedures,  or  even  to  derive 
all  of  the  properties  of  the  particular  filters  described  here.  For  read- 
rs  interested  in  additional  information,  Hamming  (1983)  gives  a  good 


a.  A 


Figure  l.  Transfer  functions 
of  the  four  types  of  filters: 
a)  lowpass,  b)  highpass,  c)  band¬ 
pass,  and  d)  bandstop 


general  introduction  to  the  topic  of  digital  filters,  while  Stearns  (1975) 
provides  more  mathematical  details.  Otnes  and  Enochson  (1978)  offer  some 
specific  examples  of  digital  filtering. 

Two  types  of  subroutines  are  described  in  this  report.  The  first  type 
designs  the  filter  (calculates  the  filter  coefficients)  and  is  denoted  by 
the  suffix  DES  in  the  subroutine  name.  The  second  type  applies  the  filter 
(convolves  the  filter  coefficients  with  the  time  series  data  to  produce  the 
output  series)  and  is  denoted  by  the  suffix  FILT.  Subroutines  LPDES, 

HPDES  and  BPDES  are  reproduced  with  slight  changes  from  Stearns  (1975). 
Subroutine  BSDES  and  all  of  the  FILT  subroutines  were  written  by  the 
author,  using  Stearns  (1975)  as  a  reference. 

The  lowpass  filter  is  based  on  the  classical  Butterworth  analog  filter 
design  procedure  (Hamming  1983,  p  222  ff.).  The  Butterworth  filter  is 
defined  by  its  power  gain 


I  H(u> ) 


1 2  = 


1  +  (— ) 


In 


where  N  is  the  order  of  the  filter  and  <uc  is  the  cutoff  or  design  fre¬ 
quency  (Fig.  2).  H(w)  is  the  transfer  function  of  the  filter.  This  filter 
has  maximal  smoothness  in  the  pass  and  the  reject  band  and  is  tangent  at 
the  origin  and  at  infinity.  By  specifying  the  amount  of  ripple  allowable 
in  the  pass  and  reject  bands,  the  order  and  the  location  of  the  poles 
needed  to  implement  the  filter  can  be  determined.  The  bilinear  transforma¬ 
tion 


(where  i  is  /-I  and  z  is  the  new  complex  variable)  is  then  used  to 
determine  the  digital  filter  parameters. 

The  other  filters  are  also  based  on  the  Butterworth  filter  design 
procedure,  except  that  a  mapping  procedure  is  used  to  convert  the  lowpass 
filter  into  the  highpass,  bandpass  or  bandstop  forms  (Stearns  1975,  Table 
12-2)  or  before  the  bilinear  transformation  is  used  to  obtain  the  digital 
form.  Details  are  given  in  Appendix  A  for  the  bandstop  filter. 

To  use  the  digital  filters,  a  design  (DES)  subroutine  is  first  called 
with  the  cutoff  frequency  (or  frequencies),  data  sample  interval  and  filter 
order  specified.  Then  a  filtering  (FILT)  subroutine  is  called  to  applv  the 
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Figure  2.  Transfer  function  of 
a  Butterworth  lowpass  filter. 

filter.  This  subroutine  convolves  the  filter  coefficients  with  the  input 
data  series  to  produce  the  filtered  output  data  series.  The  convolution  is 
then  repeated  with  the  data  series  in  reverse  order.  This- second  filter¬ 
ing  removes  any  time  delays  produced  by  the  first  filtering  process,  and 
the  net  result  is  that  the  filter  act;  as  a  zero-phase  filter  with  a  trans¬ 
fer  function  that  is  the  square  of  the  transfer  function  of  the  filter 
specified  in  the  design  stage. 

The  subroutine  parameters  are  listed  in  the  next  section.  Complete 
listings  are  given  in  Appendix  B  and  are  available  to  users  on  CRREL’s 
PRIME  9750  computer  system.  The  FILT  subroutines  require  subroutine  REVERS 
which  was  written  by  Robinson  (1978).  This  subroutine  is  also  listed  in 
Appendix  B. 

DESCRIPTION  OF  SUBROUTINES 
Filter  design  subroutines 
Lowpass  filter 

Calling  procedure:  CALL  LPDES(FC,DT,NS,A,B,C). 

Input  variables: 

FC  =  cutoff  (half-power  or  -3  dB)  frequency  in  hertz 
DT  =  time  interval  between  samples  in  seconds 
NS  =  number  of  filter  sections 
(2*NS  =  order  of  the  filter). 

Output  variables: 

A ( X ) , B ( K )  ,  C ( K  )  =  filter  coefficients  for  K  =  1,2, ...NS;  double 
precision  real  arravs. 


Highpass  filter 


Calling  procedure:  CALL  HPDES ( FC , DT ,NS , A , B ,C) . 

Input  and  output  variables:  Same  as  for  lowpass  filter  (above). 

Bandpass  filter 

Calling  procedure:  CALL  BPDES(F1 ,F2 ,DT ,NS ,A , B ,C ,D ,E ) . 

Input  variables: 

FI  =  lower  cutoff  frequency  in  hertz 
F2  =  higher  cutoff  frequency  in  hertz. 

Output  variables: 

A(K) ,  .  . . ,E(K)  =  filter  coefficients  for  K=1,2,...NS;  double  precision 
real  arrays. 

All  others  same  as  for  lowpass  filter  (above). 

Bandstop  filter 

Calling  procedure:  CALL  BSDES (FI  ,F2 ,DT ,NS , A ,B ,C ,D ,E ,F ,G) . 

Output  variables:  A(K) , . . . ,G(K)  =  filter  coefficients  for  K-1,2,...NS; 

double  precision  real  arrays. 

All  others  same  as  for  bandpass  filter  (above). 

Filtering  subroutines 

Lowpass  filter 

Calling  procedure:  CALL  LPFILT (X ,LX ,NS , A , B , C  )  . 

Input  variables: 

X  =  real  input  data  array 

LX  =  integer  number  of  entries  in  input  data  array 
NS  =  number  of  filter  sections 

A(K)  ,B(K)  ,C(K)  =  filter  coefficients  for  K=1,2,...NX;  double  precision 
real  arravs. 

Output  variables: 

X  =  real  output  data  arrav  containing  the  filtered  data. 

(Note  that  the  input  data  will  be  lost.) 

Subroutine  called:  REVERS(X,LX)  to  reverse  the  order  of  a  real  arrav 

(Robinson  1978). 

Highpass,  bandpass  and  bandstop  filters 
Calling  procedures:  CALL  HPFILT(X , LX ,NS , A  ,  B ,C ) 

CALL  RPFILT(X,LX  ,jNS,A,B,C,D,E) 

CALL  BSFILKX.LX.NS.A.B.C.D.E.F.G). 


The  number  of  filter  coefficients  changes  for  the  bandpass  and  bandstop 
filter,  otherwise  all  parameters  are  the  same  as  for  the  lowpass  filter. 

USER  INSTRUCTIONS 

To  filter  a  time  series,  the  desired  number  of  filter  sections  NS  and 
the  cutoff  frequency  (or  frequencies)  should  be  specified  in  a  call  to  sub¬ 
routine  XXDES,  where  XX  is  replaced  by  the  type  of  filter  desired  (LP,HP, 

RP  or  RS ) .  The  time  interval  DT  between  samples  is  fixed  bv  the  data  time 
series.  Frequencies  should  be  specified  in  the  reciprocal  of  the  units 
used  to  specify  DT  (e.g.,  hertz  if  DT  is  in  seconds,  cycles/day  if  DT  is  in 
davs,  etc).  The  filter  coefficients  are  returned  in  double  precision 
arravs.  These  coefficients  are  then  used  in  subroutine  XXFILT  to  carry  out 
the  actual  filtering  of  the  data  time  series. 

FILTER  PROPERTIES 

Figures  3  and  4  illustrate  the  performance  characteristics  of  some 
lowpass  filters  obtained  using  subroutines  LPDES  and  LPFILT.  The  perform¬ 
ance  characteristics  of  highpass,  bandpass  and  bandstop  filters  are  dis¬ 
played  by  the  figures  in  Appendix  C.  The  data  sampling  interval  was  0.001 
seconds,  and  in  each  case  the  input  data  series  consisted  of  511  zeroes 
with  a  single  unit-valued  spike  in  the  250th  position.  Thus,  the  filter 
output  represents  the  impulse  response  of  the  filter.  The  Fourier  trans¬ 
form  of  the  impulse  response  of  the  filter  was  then  calculated  and  the 
amplitude  spectrum,  which  is  equivalent  to  the  transfer  function  of  the 


a.  Linear  display  of  b.  Logarithmic  display 

the  transfer  function.  of  the  transfer  function. 

Figure  3.  Lowpass  filter  performance  characteristics.  One  stage  filters 
with  FC  =  400,  200,  100,  50  and  5  Hz. 
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Figure  3  (coat’d)-  Low- 
pass  filter  performance 
characteristics.  One 
stage  filters  with  FC  = 

400 ,  200,  100,  50  and  5  Hz 


filter,  was  found  and  plotted  on  a  normalized  scale.  A  logarithmic  plot  of 
the  transfer  function  is  also  given  in  the  figures. 

The  figures  show  that  the  impulse  response  of  the  filter  broadens  as 
the  passband  of  the  filter  decreases.  This  property  is  true  for  all  filt¬ 
ers,  and  the  user  will  sometimes  have  to  compromise  between  the  narrowness 
of  the  frequency  spectrum  and  the  duration  of  the  signal  in  the  time  do¬ 
main.  For  example,  Figure  3  shows  that  for  a  5-Hz  bandwidth,  the  single 
spike  that  was  originally  0.001  seconds  in  duration  is  now  0.25  seconds 
wide . 

A  comparison  of  the  plots  in  Figure  4  shows  that  as  the  number  of 
filter  sections  increases,  the  slope  of  the  transfer  function  between  the 
passband  and  the  stopband  steepens.  The  filter  performance  improves  at  the 


.  V  *ef*d  Data 


g&lllg 

Wpfly* 

jgll 

Dec-bei  Spectrum 

,  I 


1  itereJ  Dot  a 


-i  0 

4  -40 


Frequency  (Hz) 


Decibel  Spectrum 


Lowpass  Filter  FC=240  Hz  ^ 
Impulse  Response 


100  200  250 

r requeue  y  l^z 


T i me  ( s  ) 
Decibel  Spectrum 


2>0  25o 

Frequency  (Hz ) 


ire  ' .  ''x ample  of  lowpass  filtering  applied  to  real  data:  a)  unfil- 

•  ■  !  dat  i  and  decibel  amplitude  spectra,  h)  lowpass  filter  impulse 
p  vise  n|  de.-iV'l  transfer  function,  and  c )  filtered  data  and  decibel 


I 


Ur 


expense  of  increased  computational  time  and  larger  side  lobes  in  the 
impulse  response. 

Figure  5  shows  the  results  of  application  of  a  lowpass  filter  to 
actual  data.  In  digitizing  some  seismic  data,  it  was  noticed  that  the 
digitizing  process  introduced  a  strong  signal  near  the  Nyquist  frequencv  of 
250  Hz  (see  the  increase  in  decibel  spectrum  near  this  frequency  in  the 
unfiltered  data).  This  "iitter"  was  removed  by  deriving  and  applying  a 
lowpass  filter  with  a  cutoff  frequencv  of  240  Hz,  with  the  results  shown. 


SUMMARY 

FORTRAN  subroutines  that  design  and  apply  zero-phase  frequencv  filters 
to  digital  data  have  been  described.  User  instructions  have  been  included 
and  the  general  properties  of  these  filters  have  been  pointed  out. 

Complete  listings  of  the  subroutines  are  given  in  Appendix  B. 
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APPENDIX  A:  DERIVATION  OF  THE  BANDSTOP  FILTER 


The  bandstop  filter  was  derived  by  applying  a  frequency  transformation 
that  maps  a  Butterworth  lowpass  filter  transfer  function  into  a  bandstop 
shape.  The  procedure  shown  by  Stearns  (1975,  p.  201  ff.)  is  in  two  steps: 

1.  Transform  an  analog  lowpass  transfer  function  H^p(s)  into  an 
analog  bandstop  transfer  function  Hgg(S). 

2.  Apply  the  bilinear  transformat  ion 


S 


z-1 
z+ 1 


to  convert  the  analog  transfer  function  to  a  digital  transfer  function 
Hgg(z).  The  Butterworth  lowpass  analog  filter  transfer  function  is 

hlp(s)  ”  (S-S^'d-Sj) 

where  S  =  iu>  =  i 2 tt  f .  The  mapping  that  converts  this  transfer  function  into 
bandstop  form  is 

Sw2 

S  -*•  - - - 

c2  . 

S  +  uj  ^  u>2 

(see  Fig.  A]).  Doing  the  mapping  gives 


hbs(s) 


u)2(S;+S2) 

s4  -  c  _ s3  + 

Sj  s2 


+  2  S2  +  _ 

u)2  ( S  j  +  S2  ) 

[ -  +  2  a)] it>2  )  -  u)jO)2S+u^w^ 

^  1  s2  blb2 


- 1  m  r —  . 

A  .  A  A  /■  .  A 

Figure  A1 .  The  transfor¬ 
mation  S  *  S  ,\2/S2  +  -'2 

converts  the  lowpass  trans¬ 
fer  function  (a)  into  a  band 
stop  transfer  function  (b) 
(after  Stearns  1975,  Table 


If  we  select  0)1  and  0)2  to  be  the  half-power  or  -3  dB  points,  then  the 
poles  can  be  written  in  exponential  form  as 


S,  =  R  i6 
1  e 


So  =  r  -ie 

*  e 


where 


R  —  o)^  —  o)2  —  o)^ 

8  =  IT  C2n  +  N  -  1)/2N  ,  n  =  1,2,3, ...2N 

and  N  is  the  order  of  the  filter.  These  values  are  determined  by  the 
specification  of  the  Rutterworth  filter  during  the  design  phase.  Then 

Sj  +  S2  -  2  R  cos8 

S,  So  =  R2 


0)2  (  S,  +  So  ) 


2  R  cos0 


-g- V1  ■  +  2  u>iu)2  =  tof  + 

bl  b2 

Substituting  these  values  into  the  bandstop  transfer  function  gives 


"bs(s)  ■ 


4  ,  2  2 
S  +  2  u>i  u>2  S  +  coju)2 


S  -  2  R  c os 0  S  +  (u)*j  +  o)^ )  S^  -  2  R  cos0  o)jo)2  S  +  0)10)2 

Applying  the  bilinear  transformation  converts  this  analog  transfer  function 
into  a  digital  transfer  function: 


2Rcos0  (— rj  +  (W5+0)})  (— r)  -  2Rcos0  id1o)2  (— r)  +  u^o)$ 


(z-1)4  +  2o)  1 0)2  (  z- 1 ) 2  (  z+ 1  )2  +  o>^o)£(z+l) 


s:  i 


BS(z)  =  {[(l+u^)2]  (z4+l) 

+  [  -4(  1+0)^012 )( l-u)i0)2  )  ]  (z3+z2) 

r  2  2  -i  2 1 

+  [6  +  6o)ju>2  -  4u>iU)2jz  } 

22 

/ { [  1  +  0)^012  +  2  W2  ~  2(l+u>iu>2)  R  cos9]  z 

2  2  i 

+  [-4+4  0)i0)2  +  4(1-<i)j0)2)  R  cos  0]  zi 

r  22  i  2 

+  [6+6  ujjW2  “  4  uii oi2 J  z 

2  2 

+  [-4+4  0)^2  ~  4  (l-o)ju)2)  R  cos  0]  z 
2  2 

+  [  l+a>i(i)2  +  2co i CO2  +  2  (l+0)ju)2)  R  cos©]} 

A  digital  bandstop  filter  with  the  above  transfer  function  was  implemented 
by  the  FORTRAN  subroutines  BSDES  and  BSFILT. 


APPENDIX  B:  LISTING  OF  FORTRAN  SUBROUTINES 


The  following  suhriiut i nes  for  digital  filtering  are  listed  in  this 
appendix:  LPDES ,  HPDES ,  BPDES ,  RSDES ,  LPFILT ,  HPFILT ,  BPFILT ,  BSFILT, 

REVERS. 
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SUBROUTINE  '_?DESirC,DT  NS  At3  •:> 
i  PASS  307”EawCR-H  DiSIViL  F  t._  fZR  DESIGN  SBR 


REFERENCE.*  S.D.  S "EARNS  (iiS/5)  DIGITAL  SIGNAL  ANALYSIS,  P.  2 S3. 


INPUTS: 

FC  =  CUTOFF  FREQUENCY  1-3  09)  IN  HZ 
DT  =  SAMPLING  INTERVAL  IN  SECONDS 
NS  =  NUM9ER  OF  f;lT£R  SECTIONS 
OUTPUTS: 

A i K  j , 9 ! K  ) ,  AND  C!K),  FILTER  COEFFICIENTS  FOR  <-1  TO  NS 
Th£  DIoITAl  FIlTER  HAS  tS  SECTIONS  IN  CASCADE.  THE  KTH 
SECTION  HAS  THE  TRANSFER  FuNCTION 


H(Z)  = 


At  K  >*( Z**2  +  2*Z  +1  ) 
Zx*2  *  B<K)xZ  *  C(K) 


IF  F(M)  AND  G(M)  ARE  THE  INPUT  AND  THE  OUTPUT  OF  THE  KTH 
SECTION  AT  TIME  M*T,  THEN 


GtM)  =  AtO*(F(M)  +  2xF(  M- 1  )  +  F(M-2)) 
-9( <  )*G: M- !  )  ~Ct  K ) *G( M-2 ) 


REAL*8  A< i  ),3( 1), C( t  ) 

PI  =  3. 1 4 1 5826536 
ALPHA  =  FC*PI*T 
WCP  =  TAN (ALPHA) 

WCP2  =  WCPxvkCP 
DO  10  K  =  I,  NS 

CS  =  C0S(^L0AT(2*(K+NS)-n*PI/FL0AT<4*NS)  ) 
X  =  1 ./<  t. *WCP2~2*VCP*CS > 

A( X  )  =  WCP2*X 
BtX  )  =  2. *<  WCP2-1 .  !*X 
CiK)  a  ( 1 . +WCP2+2.*WCP*CS)*X 
CONTINUE 
RETURN 
END 


******xxxx**xxx*****xxxxx**x**x**x***x********x****x** 

SUBROUTINE  HPDES ( FC, DT, NS. A,B,C) 

HIGH  PASS  BUTTERVORTH  &IG1TAL  £ ILTER  DESIGN  SBR 

REFERENCE;  S.D.  STEARNS  i!875)  DIGITAL  SIGNAL  ANALYSIS,  P.  270. 

INPUTS; 

FC  =  CUTOFF  FREQUENCY  (-3  DB)  IN  HZ 
DT  =  SAMPLING  INTERVAL  IN  SECONDS 
NS  =  NUMBER  OF  FILTER  SECTIONS 
OUTPUTS: 

At  K  J , Bt  K  J ,  AND  C ( K  J ,  FILTER  COEFFICIENTS  FOR  K=1  TO  NS 
THE  DIGITAL  FILTER  HAS  N§  SECTIONS  IN  CASCADE.  THE  KTH 
SECTION  HAS  THE  TRANSFER  FUNCTION 


H(Z>  = 


At  K )*t  Zx*2  -  2xZ  +1  ) 
Zx*2  *  8t  K  )*Z  +  C(K) 


IF  FtM)  AND  G(M)  ARE  THE  INPUT  AND  THE  OUTPUT  OF  THE  KTH 
SECTION  AT  TIME  M*T,  THEN 


GlM)  =  At  K  )*( Ft  N  )  -  2*F(M-1  )  +  F(M-2)) 
-3(K)*G<M- t  )  -C<K)*G1M-2) 


10 


REAL*8  At  1  >,3( !  ),C(1) 

P  1  =  3.  M1S92S536 
*CP=SIN(FCxPI*T)/COS(FCxP;*T) 

DO  10  K=1,NS 

CS=CCS( £lCAT ( 2*<K*NS  >-!  )*P I /FLOAT ( 4*NS) ) 
At  K  )  = ! ./! 1 . *WCP*wCP-2 . *WCP*CS) 

B<K  )=2.*(WCP*WCP-1 ) *A(  K ) 

C( K  )  =  ( 1 . -*C?xWCP>2. *WCP*CS)*A(K ) 

CONTINUE 

RETURN 

END 


***xx* *x**x* x *** *xx **x ********** ********* 


nnonnnnnnnnoonnonn(  »r  ion 


C  a**************************************** 

SL3R0uT;nE  3PDESI  r '■  F2,DT,NS,A3,C,D,E) 

SANDBAGS  9UT 'ESVCRTH  DIGITAL  ~ IL ;>£s!GN  S3R 
REFERENCE:  S.O.  STEARNS  (IS75)  DIGiTAL  SIGNAL  ANALYSIS,  P. 
I NPLT3: 

■:,r2  =  r.jTCFF  FREGL'ENC  IES  (-3  DB  >  IN  HZ 
DT  =  SAMPLING  interval  IN  SECONDS 
NS  =  NJH3ER  OF  FIL.ter  SECTIONS 
OUTPUTS: 

A-.  K>  THROUGH  ECO,  filter  COEFFICIENTS  FCR  K  =  1  TO  NS 
T'-£  DIGITAL  filter  -as  ns  sections  in  cascade,  the  <th 
SECTION  HAS  the  transfer  function 


A(<  !*(Z**4  -  2*Zt*2  *t ) 

HiZ)  =  - 

Z**4  *  BlX)*Z**3  ♦  Ci  K  )*Z**2  *  D<K)*Z  ♦  El K ) 

IF  F(H)  AND  3(M)  ARE  THE  INPUT  AND  THE  OUTPUT  OF  THE  KTH 
SECTION  AT  TIME  M*T,  THEN 

G(M)  =»  A(K)*lF(M>  -  2*F  ( M-2  )  ♦  F(M-4)> 

-3( K )*Gl M-1  )  -C(  K  1 *G<  H-2 ) 

-Di K ) *G<  M-3 )  -E( K )*Gl M-4 ) 

REAL*8  A( 1  ),8< 1  ),C( 1  ),D( 1 >,EC I) 

PI=3. 141592^536 

W 1 =S INI F 1 *PI*T ) /COSl FI *PI*T ) 

W2=SINlF2*PUT)/C0S(F2*PUT) 

WC=W2-W1 

G=WC*WC+2.*W1*W2 
S=W1 *WI *W2*W2 
DO  ! 0  K= 1 , NS 

CS=C0SlrL0AT(2*lK+NS)-1  )*PI/FLOAT(  4*NS  )  ) 

P=-2 . *WC*CS 
R=P*V  UW2 

X=1  .  *P  +  Q  +  R«-S 
Al K )=WC*WC/X 

8(  K  )  =  (  -  4.  -  Z.*P  *■  2.  *R  *  4.  *S)/X 
C(  K  )  =  (  6.  -  2.  *Q  *  6.  *S  >/X 
D( K  )  =  (  -  4.  v  2.*P  -  2. *R  ♦  4 . *S  )/X 
E<  K  )  =  <  1  .  -  P  *  0  -  R  *  S)/X 
10  CONTINUE 
RETURN 
END 


271 


ooonnooononoonnnoonnnonnno  n 


************************ ****** ******** 

SUBROUTINE  SSDESlFI ,F2,DT,NS,  A, 3, C^D , E, F , G > 
BANDS  TC?  3UTTERWCRTH  DIG!  tAL  FlC^R  D^So^  $3R 
THIS  S3R  DERIVED  AND  DESIGNED  BASED  DN  THE  BOOK 
S.D.  STEARNS,  1975,  D I G I TAL  SIGNAL  ANALYSIS, 
INPUTS: 

'  ,F2  =  CUTOFF  FREQUENCIES  (-3  DB)  IN  HZ 
=  SAMPLING  INTERVAL  IN  SECONDS 
=  NUMBER  OF  FILTER  SECTIONS 


HAYDEN  800K  CO. 


c  < 

dt’ 

NS 

OUTPUTS: 

AIK)  THROUGH  EIK),  filter  coefficients  for  K=l  TO  NS 


THE  DIGITAL  FILTER  has  ns  sections  in  cascade. 
HAS  THE  TRANSFER  FUNCTION 


THE  KTH  SECTION 


HIZ)  = 


A( K  )*( Z**4  +  n  ♦  B( K )*( Z**3  *  Z)  +  C<K)*Z**2 


Z**4  *  D(K)*Z**3  *  El K ) *Z**2  +  F(K)*Z  +  GlK> 


IF  F( M )  AND  G(M)  ARE  THE  INPUT  AND  OUTPUT  OF  THE  KTH  SECTION 
AT  TIME  M*T,  THEN 


GIM)  =  Al K )*( F( M )  +  F I M-4 ) ) 

♦Bl K  )*( F I M- 1  )  +  F ! M-3 )  )  +  Ct K )*F( M-2 ) 

-D(K)*G(M-1)  -  EIK) *G( M-2 )  -  F(K)*G(M-3>  -  G(K)*G<M-4> 


USE  SBR  BSFILT  TO  APPLY  THE  8ANDST0P  FILTER 
DON  ALBERT  12/21/82 


1  MLOCn  I  \  I  / 

REAL*8  A( 1 > ,  B< 1 ),C( 1 ),D( 1 ),E< 1 ),F(1 ) ,  3(  1 ) 
P I  =3 . 141 5926536 


<  *  —  •  I  T  I 

W1=SIN(F1*PI*T) /COS  I  F 1 *PI*T ) 
W2=SINIF2*PI*T)/C0S(F2*PI*T) 

Q=*I1  *W2 
Q2=Q*Q 

qp=i . +q 

QM=1 . -Q 
R=W2-W1 

S=WUWt  +  V2*V2 
DO  10  K=1,NS 

CS=R*COSl FLOAT! 2*1 K+NS>-1 )*PI /FLOAT!  4*NS > ) 
X»1./(1.+Q2+S-2.*CS*QP> 

A I K  )=QP*QP*X 
B<K>=-4.*QP*QM*X 
Cl  K  )  =  (  6.  +S. *Q2-4 . *0 )  *X 
D<  K  )=4.  *1  -1  .  *QP*QMt-CS*QM) 
E<K)=I6.+S.*Q2-2.*S> 

FI K  )  -4 . *( - 1 . *QP*QM-CS*QM) 

G( K )=< I . +Q2+S+2. *CS*QP ) 

CONTINUE 
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X  IXXIlllltXX^XXXXtttl^iXtttlXtt^xXXOXltXlXlltX 

SUBROUTINE  -PFILTc X  '  X. VS. A, 3. C) 
applies  A  ..3W  PASS  =:_ter  TO  ARR^Y  X 
calls  ssr  REvERS 
BASED  ON  INFO  CN  STEARNS  (1975) 
don  assert  8/23/83 


REAL *8  A( 1  ), B< 1  ),c; 1 > 
.REAL*4  X (  LX  ) 

RE AL*4  Y(  >1,3) 

RANDLES  FIl'EF.S  UP 


.  _  _  _  TO  NS=i0 

APPLIES  THE  FILTER  'VICE,  ONCE  IN  EACH  DIRECTION. 
I F:_AG=3 
CCNT INLE 
I~LA3=IFLAGH 
I  F( I FLAG. EG. 3)  RETURN 
NS  1  =NSM 
DO  10  I  =  1 , NS  t 
DO  I  0  J= 1 ,3 
Y( I . J)=0.0 
CONTINUE 
DO  40  Hal. LX 
Y ( 1,3)  =X( M ) 

DO  £0  N= I , NS 

Y  ( N+1 ,3)=A(N)*(Y(N.3)+2.*Y(N,2)+Y(N,  1  )) 
-B(N)*Y(N+!  ,2)  -  C(N)*Y<NH  ,  1  ) 

CONTINUE 
DO  30  N= 1 . NS  1 
DO  30  mN=1 , 2 
Y(  N. MM ) =Y<  N, MM+ 1  ) 

CONTINUE 
Xl  M  )  =  Y<  NS !  ,3) 

CONTINUE 

CALL  REVERS( X, LX ) 

GOTO  5 
END 


I 


************************* ********************** 

SUBROUTINE  HPFILHX.LX.NS,  A.B.C) 
applies  a  high  PASS  FlLTEfc  T0  ^rray  X 
calls  sbr  reyers 

BASED  ON  INFO  IN  STEARNS  (I97S) 

DON  ALBERT  6/30/83 


REAU*8  A< 1  ),B( I  ),C< I ) 

REAL*4  X ( LX ) 

REAL*4  Y ( I  1,3) 

HANDLES  filter^  UP  TO  NS- 1 0  „  ,„kl 

APPLIES  the  FILTER  twice,  once  in  each  direction, 
IFLAG=0 

continue 

Tflag=;flag+i 

IF  (  IF^.AG.EQ.3)  RETURN 
NS  I  =NSM 
DO  '0  1  =  1, NS  I 
00  10  0=1,3 
Y ( I, J)=0.0 
CONTINUE 
DO  40  M= 1 , LX 
Y ( 1 , 3)*X(M) 

DO  20  N= I , NS 


20 


1 


Y<NH  ,  3)=A(N)*(  Y<N,3)-2.  *Y<  N,  2  )«-Y<  N,  I  >  > 
-B(N)*Y(N+I,2)  -  C(N)*Y(N-H  ,  !  ) 


30 

*0 


CONTINUE 
DO  30  N= 1 , NS  1 
DO  30  MN=1 . 2 

Y!N.*M>=Y<N,MM*1 > 
CONTINUE 
X<  M  )=Y<  NS1 ,3) 

continue 

CALu  SEVERS! X, LX) 

GOTO  5 
END 


SO  the  PHASE  IS  2ERC 


SO  THE'  PHASE  IS  ZERO 


SUBROUTINE  5PFILT <  X  LX, NS, A, 3, C, D,  E ) 

APPLIES  A  3ANDPASS  FILTER  TO  A^rAy  X 

CALLS  SBR  REVERS 

BASED  ON  INFO  !N  STEARNS  1)875) 

DON  ALBERT  1 2/21/32 

REAL  *8  A(  I  ).  BU  ),C(  t  ),D(I  ).EU  > 

REAL*4  X',  LX) 

R£Al*4  YU  1  5  ) 

HANDLES  FILTERS  JP  TC  NS=13 

APPLIES  'HE  FILTER  T*IC£,  ONCE  IN  EACH  DIRECTION,  SO 
:PuAG=0 
CONTINUE 
I^k-AO  =  ;Fi_AG+ 1 
IFC I FLAG. £Q. 3 )  RETURN 
NSI =NS+1 
DO  13  1=1, NS! 

DC  10  U=  1  ,4 
Y(  l.  J)=0.6 
CONTINUE 
DC  40  M=1 , LX 
Y ( 1 ,5)=X(h) 

DO  ^0  N=) , NS 

Y(N+1 , 5)=A(N)*<  Y<N,5)-2.*YCN, 3>+Y(N, t ) ) 

1  ~B( N ) *T ( N* 1 , 4)-C(N)*YlN+1  ,  3) 

2  -Di  N )*Y( N+ 1 , 2 )-£t  N  )*Y( N+ 1 ,  1  ) 

CONTINUE 

DO  30  N=l, NSI 
DO  30  M^i=l  ,  4 
YIN, MM )=Y<  N, MM+ 1  ) 

CONTINUE 
X(M)=Y(NS! ,5) 

CONTINUE 

CALL  REVERSl X, LX ) 

GOTO  S 
END 


HE  PHASE  IS  ZERO. 


********************************************************* 

SUBROUTINE  8SFILT  ^  X, LX, NS, A, B, C,  D,  E,  F,  G  ) 

APPLIES  A  BANDSTOP  FrLfER  TO  AI^RAY  X 

CALLS  SSR  REVERS 

BASED  ON  INFO  IN  STEARNS  0  875) 

CON  AL3ERT  12/21/82 

REAL*8  A( 1  ),B( t  ),C( I  ),D< 1  1 , E( 1 ) ,  F( 1  ),G( 1 ) 

REAL*4  X<LX) 

REAL*4  Ylll.Sl 

HANDLES  FILTER^  UP  TO  NS*  1 0  te. 

APPLIES  THE  FILTER  tvicE,  ONCE  IN  EACH  DIRECTION,  SO  THE  PHASE  IS  ZERO. 
JFLAG=9 
CONTINUE 
IFLAG=IFLAG*1 
IF( IFlAG. EQ. 3 )  RETURN 
NS1=NS+1 
DO  10  1=1, NSI 
DO  13  J=1,5 
Y (  I,  J 1=0. ^ 

CONTINUE 
DO  40  M=  t , LX 
Y (  1 , 5  )  =X(  M  ) 

°°YfN+,1<rs')=A(N)*(  Y<N,5)+Y<N,  )))  *•  B<  N  )*<  Y<  N^4  )+Y(  N,  2  >  > 

1  +C(N)4Y(N,  3) 

2  -D(N!*Y(N+1 ,4)  -E(N)*Y(N«-!,3> 

3  -F( N ) *Y ( N* 1 , 2 )  -GIN )*Y< N*1 , 1) 

CONTINUE 

DO  30  N=1 , NSI 
DC  30  Mn= I . 4 
Y . N . MM , =Y( N, MM*1  I 
CONTINUE 
XIM)=Y(NS1 .5) 

CONTINUE 

CALL  REV£RS( X, LX) 

GOTO  5 
END 

*****************  *************************** ***a**as.*i* 

SUSTCu'INE  REvERS( X, LX ) 


REVERSES  A  TIME  SERIES 
REFERENCE:  ROBINSON  < 1 5S7 ) . 
DIMENSION  Y(LX) 

L=LX'2 
DO  13  I  =  1  ,  L 
J=LX-I 
TEMF=X< !  ) 

X(  I  >=X(  >1  ) 

X( J+'  > =  EYP 


P  28. 


APPENDIX  C:  FILTER  PERFORMANCE  CHARACTERISTICS 


For  each  figure  in  this  appendix,  a  filter  was  applied  to  an  input 
time  series  consisting  of  511  zeroes  with  a  unit-valued  spike  at  the  250th 
position  to  produce  the  filter's  impulse  response.  The  data  sampling 
interval  was  0.001  seconds.  The  Fourier  transform  of  the  impulse  response 
was  then  calculated  and  used  to  determine  the  filter's  transfer  function, 
which  was  plotted  on  a  normalized  linear  and  on  a  logarithmic  scale. 


Amplitude 


a.  Linear  display  of 
the  transfer  function. 


Decibel 


b.  Logarithmic  display 
of  the  transfer  function. 


c.  Impulse  response  (actual  data  are 
symmetric  about  the  center). 


Impulse  Response 


r  . . . 

! 

1 

i _ 

100 

i 

i 

- N 

- - - 

l 


50 _ | 


] 


0 


Time 


(  s) 


1 _ I 

0  511 


Figure  Cl.  Highpass  filter  performance  characteristics.  One  stage  filter 
with  FC  =  400,  200,  100,  50  and  5  Hz. 


Impulse  Response 


0  0  511 

Time  (s  ) 


c-  Impulse  response. 

Figure  C6  •  Bandpass  filter  performance  characteristics.  One,  two  and 
four  stage  filters  with  lower  FC  =  100  Hz,  upper  FC  =  150  Hz. 


